ETC1010/5510 Tutorial 11 - Solution

Introduction to Data Analysis

Author

Patrick Li

Published

October 7, 2024

🎯 Workshop Objectives

  • Introduction to linear models
  • Interpreting linear models
  • Linear model diagnostic

🔧 Instructions

Install the necessary pacakges

🌶️ Exercise: Introduction to Linear Model

library(tidyverse)
library(broom)

Section A: Practice with Advertising Data

1. Read the data into RStudio. Remove the first column.

advert <- read_csv("data/Advertising.csv") %>% 
  select(-1)

head(advert)
# A tibble: 6 × 4
     TV radio newspaper sales
  <dbl> <dbl>     <dbl> <dbl>
1 230.   37.8      69.2  22.1
2  44.5  39.3      45.1  10.4
3  17.2  45.9      69.3   9.3
4 152.   41.3      58.5  18.5
5 181.   10.8      58.4  12.9
6   8.7  48.9      75     7.2

Before fitting a model, let’s do some EDA. Plot the data with sales on the y-axis, and the rest of the variables on the x-axis in 3 separate plots. Then combine them into one plot.

p1 <- ggplot(advert,
       aes(x = radio,
           y = sales)) +
  geom_point(alpha = 0.2)

p2 <- ggplot(advert,
       aes(x = newspaper,
           y = sales)) +
  geom_point(alpha = 0.2)

p3 <- ggplot(advert,
       aes(x = TV,
           y = sales)) +
  geom_point(alpha = 0.2)

gridExtra::grid.arrange(p1, p2, p3, nrow = 1)

2. Fit a model with just radio as the x variable.

The simple linear regression is of the form:

\[y_i=\beta_0+\beta_1x_{i}+\varepsilon_i, ~~~ i=1, \dots, n\]

Suppose we are interested in the amount of variation in sales that is explained by the variation in radio.

advert_m1 <- lm(sales ~ radio, data = advert)
advert_m1

Call:
lm(formula = sales ~ radio, data = advert)

Coefficients:
(Intercept)        radio  
     9.3116       0.2025  

The formula response ~ explanatory specifies the response variable and explanatory variable from the data. In this case, the response is sales in in thousands of units and the explanatory variable is spending on radio advertisement, in thousands of dollars.

Note that by default, the lm() function will include the intercept term.

We can extract components of the regression line using the broom package. The function tidy() will extract the estimated coefficients as a tidy table.

library(broom)
coef_advert_radio <- tidy(advert_m1)
coef_advert_radio
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    9.31     0.563      16.5  3.56e-39
2 radio          0.202    0.0204      9.92 4.35e-19

The fitted model has the following functional form:

\[\widehat{\text{sales}} = 9.311 + 0.202~\text{radio}\]

The coefficients can be interpreted as

Slope: For each additional thousand dollars spending on radio advertising, on average, the sales increases by 202 units.

Intercept: When there is no advertising spending on radio, sales increases by 9 thousand units.

We can draw our fitted linear model with geom_smooth(method = "lm")

ggplot(advert,
       aes(x = radio,
           y = sales)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "lm")

3. Diagnostics and augment

The residuals from the model are the differences between the observed values and the fitted values. The residuals and fitted values can be extracted from using augment() function as part of the broom package:

advert_m1_diagnostics <- augment(advert_m1)
advert_m1_diagnostics
# A tibble: 200 × 8
   sales radio .fitted   .resid    .hat .sigma     .cooksd .std.resid
   <dbl> <dbl>   <dbl>    <dbl>   <dbl>  <dbl>       <dbl>      <dbl>
 1  22.1  37.8   17.0    5.13   0.00982   4.27 0.00722         1.21  
 2  10.4  39.3   17.3   -6.87   0.0109    4.26 0.0143         -1.62  
 3   9.3  45.9   18.6   -9.31   0.0167    4.23 0.0409         -2.20  
 4  18.5  41.3   17.7    0.825  0.0124    4.29 0.000237        0.194 
 5  12.9  10.8   11.5    1.40   0.00854   4.28 0.000467        0.329 
 6   7.2  48.9   19.2  -12.0    0.0200    4.20 0.0822         -2.84  
 7  11.8  32.8   16.0   -4.15   0.00707   4.28 0.00339        -0.975 
 8  13.2  19.6   13.3   -0.0806 0.00531   4.29 0.000000952    -0.0189
 9   4.8   2.1    9.74  -4.94   0.0152    4.27 0.0105         -1.16  
10  10.6   2.6    9.84   0.762  0.0147    4.29 0.000241        0.180 
# ℹ 190 more rows

We can visualise the residuals directly on the scatter plot:

ggplot(advert,
       aes(x = radio,
           y = sales)) +
  geom_point(alpha = 0.2)+
  geom_smooth(method = "lm", se = FALSE) +
  # overlay fitted values
  geom_point(data = advert_m1_diagnostics, 
             aes(y = .fitted), 
             color = "blue", 
             alpha = 0.2) +
  # draw a line segment from the fitted value to observed value
  geom_segment(data = advert_m1_diagnostics, 
               aes(xend = radio, y = .fitted, yend = sales),
               color = "blue",
               alpha = 0.2)

Another way of visualising residuals is to produce a scatterplot of the residuals against the fitted values:

ggplot(advert_m1_diagnostics,
       aes(x = .fitted, y = .resid)) +
  geom_point()

Yet, another and a better way to visualise residuals is using resid_panel() which gives you a panel of diagnostic residual plots.

library(ggResidpanel)
resid_panel(advert_m1, plots = "all")

Recall the assumptions of the linear model from the lecture. What are they?

  • Are they met from the plots?
  • Are there observations with large positive residuals? What do they mean?

4. How good is your fit?

  • \(R^2\) is a common measurement of strength of linear model fit.
  • \(R^2\) tells us % variability in response explained by the model.

We can extract model fit summaries using glance()

advert_m1_summary <- glance(advert_m1)
advert_m1_summary
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.332         0.329  4.27      98.4 4.35e-19     1  -573. 1153. 1163.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

The \(R^2\) value here can be interpreted as the proportion of variability in sales explained by spending on radio advertisement. That is approximately 33% of variability in sales is explained by spending on radio advertisement. Hmm, that is not much! 😒

Section B: Multiple Linear Regression

In reality, we will not only spend all the advertising spending on just one medium. We may choose to have other media of advertisement too such as newspaper and TV.

1. Run a multiple regression model to show the effect of different advertising media on sales. Obtain a summary of the regression results.

advert_m2 <- lm(sales~ TV + newspaper + radio,
                data = advert)
summary(advert_m2)

Call:
lm(formula = sales ~ TV + newspaper + radio, data = advert)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8277 -0.8908  0.2418  1.1893  2.8292 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.938889   0.311908   9.422   <2e-16 ***
TV           0.045765   0.001395  32.809   <2e-16 ***
newspaper   -0.001037   0.005871  -0.177     0.86    
radio        0.188530   0.008611  21.893   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.686 on 196 degrees of freedom
Multiple R-squared:  0.8972,    Adjusted R-squared:  0.8956 
F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

2. Is the relationship linear?

advert_m2_aug <- augment(advert_m2)
ggplot(data = advert_m2_aug,
       aes(x = .fitted,
           y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, colour = "red")

3.Residual Plots

Residual plots can be used in order to identify non-linearity or other problems with our regression model. * Plot the residuals using resid_panel(). Take this out of the exercise, but leave in the solutions

library(ggResidpanel)
resid_panel(advert_m2, plots = "all")

  • Do you see any patterns? [ANS: Looks fine but the residuals tends to be more positive when the fitted values are at the lower or higher values.]

  • We can see that the residuals are not normal, in which the Q-Q plot, histogram and boxplot show data skewness.

  • The Cook’s Distance plot shows our data is having two extreme outliers. You may want to investigate on this.

👴 Exercise 9B: Data modelling using the gapminder data

Gapminder

  • Hans Rosling was a Swedish doctor, academic and statistician, Professor of International Health at Karolinska Institute. Sadly he passed away in 2017.
  • He developed a keen interest in health and wealth across the globe, and the relationship with other factors like agriculture, education, and energy.
  • You can play with the gapminder data using animations at https://www.gapminder.org/tools/.

R package gapminder contains subset of the data on five year intervals from 1952 to 2007.

library(gapminder)
glimpse(gapminder)
Rows: 1,704
Columns: 6
$ country   <fct> "Afghanistan", "Afghanistan", "Afghanistan", "Afghanistan", …
$ continent <fct> Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, …
$ year      <int> 1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, …
$ lifeExp   <dbl> 28.801, 30.332, 31.997, 34.020, 36.088, 38.438, 39.854, 40.8…
$ pop       <int> 8425333, 9240934, 10267083, 11537966, 13079460, 14880372, 12…
$ gdpPercap <dbl> 779.4453, 820.8530, 853.1007, 836.1971, 739.9811, 786.1134, …

Section A: Australia

1. How has life expectancy changed over the years, for each country?

Use a line plot to answer this question and provide your explanation here.

  • There generally appears to be an increase in life expectancy in the different countries.
  • A number of countries have big dips from the 1970s through 1990s.
  • A cluster of countries starts off with low life expectancy but ends up closer to the highest life expectancy by the end of the period.

2. How has life expectancy in Australia changed over the years.

Use a line plot to answer this question and provide your explanation here.

oz <- gapminder %>% 
  filter(country == "Australia")
oz
# A tibble: 12 × 6
   country   continent  year lifeExp      pop gdpPercap
   <fct>     <fct>     <int>   <dbl>    <int>     <dbl>
 1 Australia Oceania    1952    69.1  8691212    10040.
 2 Australia Oceania    1957    70.3  9712569    10950.
 3 Australia Oceania    1962    70.9 10794968    12217.
 4 Australia Oceania    1967    71.1 11872264    14526.
 5 Australia Oceania    1972    71.9 13177000    16789.
 6 Australia Oceania    1977    73.5 14074100    18334.
 7 Australia Oceania    1982    74.7 15184200    19477.
 8 Australia Oceania    1987    76.3 16257249    21889.
 9 Australia Oceania    1992    77.6 17481977    23425.
10 Australia Oceania    1997    78.8 18565243    26998.
11 Australia Oceania    2002    80.4 19546792    30688.
12 Australia Oceania    2007    81.2 20434176    34435.

3. Fit the linear model just for Australia: lifeExp ~ year

ggplot(data = oz, 
       aes(x = year, 
           y = lifeExp)) + 
  geom_line() 

4. Fit a linear model: lifeExp ~ year.

oz_lm <- lm(lifeExp ~ year, data = oz)

oz_lm

Call:
lm(formula = lifeExp ~ year, data = oz)

Coefficients:
(Intercept)         year  
  -376.1163       0.2277  
summary(oz_lm)

Call:
lm(formula = lifeExp ~ year, data = oz)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0250 -0.5201  0.1162  0.3781  0.7909 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -376.11630   20.54716  -18.30 5.09e-09 ***
year           0.22772    0.01038   21.94 8.67e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6206 on 10 degrees of freedom
Multiple R-squared:  0.9796,    Adjusted R-squared:  0.9776 
F-statistic: 481.3 on 1 and 10 DF,  p-value: 8.667e-10

5. Tidy the model.

tidy(oz_lm)
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept) -376.      20.5        -18.3 5.09e- 9
2 year           0.228    0.0104      21.9 8.67e-10

\[\widehat{lifeExp} = -376.1163 + 0.2277~year\]

6. Center year

Let us treat 1950 as the first year, so for model fitting, we are going to shift the years to begin in 1950. This makes interpretation easier.

gap <- gapminder %>% 
  mutate(year1950 = year - 1950)

oz <- gap %>% 
  filter(country == "Australia")

7. Model for centered year

oz_lm <- lm(lifeExp ~ year1950, data = oz)

oz_lm

Call:
lm(formula = lifeExp ~ year1950, data = oz)

Coefficients:
(Intercept)     year1950  
    67.9451       0.2277  

8. Tidy the model centered year

tidy(oz_lm)
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   67.9      0.355      192.  3.70e-19
2 year1950       0.228    0.0104      21.9 8.67e-10

\[\widehat{lifeExp} = 67.9 + 0.2277~year1950 \]

9. Extract residuals and fitted values using augment()

oz_aug <- augment(oz_lm, oz)
oz_aug
# A tibble: 12 × 13
   country   continent  year lifeExp      pop gdpPercap year1950 .fitted  .resid
   <fct>     <fct>     <int>   <dbl>    <int>     <dbl>    <dbl>   <dbl>   <dbl>
 1 Australia Oceania    1952    69.1  8691212    10040.        2    68.4  0.719 
 2 Australia Oceania    1957    70.3  9712569    10950.        7    69.5  0.791 
 3 Australia Oceania    1962    70.9 10794968    12217.       12    70.7  0.252 
 4 Australia Oceania    1967    71.1 11872264    14526.       17    71.8 -0.716 
 5 Australia Oceania    1972    71.9 13177000    16789.       22    73.0 -1.02  
 6 Australia Oceania    1977    73.5 14074100    18334.       27    74.1 -0.604 
 7 Australia Oceania    1982    74.7 15184200    19477.       32    75.2 -0.492 
 8 Australia Oceania    1987    76.3 16257249    21889.       37    76.4 -0.0508
 9 Australia Oceania    1992    77.6 17481977    23425.       42    77.5  0.0505
10 Australia Oceania    1997    78.8 18565243    26998.       47    78.6  0.182 
11 Australia Oceania    2002    80.4 19546792    30688.       52    79.8  0.583 
12 Australia Oceania    2007    81.2 20434176    34435.       57    80.9  0.310 
# ℹ 4 more variables: .hat <dbl>, .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>
oz_aug$.fitted
 [1] 68.40051 69.53913 70.67775 71.81637 72.95499 74.09361 75.23223 76.37084
 [9] 77.50946 78.64808 79.78670 80.92532
oz_aug$.resid
 [1]  0.71948718  0.79086830  0.25224942 -0.71636946 -1.02498834 -0.60360723
 [7] -0.49222611 -0.05084499  0.05053613  0.18191725  0.58329837  0.30967949

10. Plot the year and lifeExp with points and add in the line.

ggplot(data = oz_aug, 
       aes(x = year, 
           y = .fitted)) + 
  geom_line(colour = "blue") + 
  geom_point(aes(x = year,
                 y = lifeExp))

11. Plot residuals against fitted values to reveal problems with the fit and interpret it.

ggplot(oz_aug,
       aes(x = .fitted, y = .resid)) +
  ylim(c(-5,5)) +
  geom_point()

# Another way to look at Residuals: .std.resid with year:
ggplot(data = oz_aug, 
             aes(x = year, 
                 y = .std.resid)) +
  ylim(c(-5,5)) + 
  geom_hline(yintercept = 0,
             colour = "white", 
             size = 2)  +
  geom_line() 

  • Life expectancy has increased 2.3 years every decade, on average.
  • There was a slow period from 1960 through to 1972, probably related to mortality during the Vietnam war.

Note: We use standardized residuals here is to make the residuals unitless. Therefore, it is easier for us to spot any unusual observations, which are those above the line of 2 and below -2. Remember the rules of 95% confidence interval?

Section B: Other countries??

1. Can we fit a model for New Zealand?

nz <- gap %>% 
  filter(country == "New Zealand")

nz_lm <- lm(lifeExp ~ year1950, data = nz)
nz_lm

Call:
lm(formula = lifeExp ~ year1950, data = nz)

Coefficients:
(Intercept)     year1950  
    68.3013       0.1928  

2. Can we fit a model for Japan?

japan <- gap %>%
  filter(country == "Japan")

japan_lm <- lm(lifeExp ~ year1950, data = japan)
japan_lm

Call:
lm(formula = lifeExp ~ year1950, data = japan)

Coefficients:
(Intercept)     year1950  
    64.4162       0.3529  

3. Can we fit a model for Italy?

italy <- gap %>%
  filter(country == "Italy")

italy_lm <- lm(lifeExp ~ year1950, data = italy)
italy_lm

Call:
lm(formula = lifeExp ~ year1950, data = italy)

Coefficients:
(Intercept)     year1950  
    66.0574       0.2697  

4. Is there a better way? What if we wanted to fit a model for ALL countries?

YES, let’s begin. 🙂

  1. Nest country level data (one row = one country)
by_country <- gap %>% 
  select(country, year1950, lifeExp, continent) %>%
  group_by(country, continent) %>% 
  nest()

by_country
# A tibble: 142 × 3
# Groups:   country, continent [142]
   country     continent data             
   <fct>       <fct>     <list>           
 1 Afghanistan Asia      <tibble [12 × 2]>
 2 Albania     Europe    <tibble [12 × 2]>
 3 Algeria     Africa    <tibble [12 × 2]>
 4 Angola      Africa    <tibble [12 × 2]>
 5 Argentina   Americas  <tibble [12 × 2]>
 6 Australia   Oceania   <tibble [12 × 2]>
 7 Austria     Europe    <tibble [12 × 2]>
 8 Bahrain     Asia      <tibble [12 × 2]>
 9 Bangladesh  Asia      <tibble [12 × 2]>
10 Belgium     Europe    <tibble [12 × 2]>
# ℹ 132 more rows
  1. What is in the data column?
by_country$data[[1]]
# A tibble: 12 × 2
   year1950 lifeExp
      <dbl>   <dbl>
 1        2    28.8
 2        7    30.3
 3       12    32.0
 4       17    34.0
 5       22    36.1
 6       27    38.4
 7       32    39.9
 8       37    40.8
 9       42    41.7
10       47    41.8
11       52    42.1
12       57    43.8

It’s a tibble! 🤪

Let’s fit a linear model to each one of them.

lm_afganistan <- lm(lifeExp ~ year1950, 
                    data = by_country$data[[1]])
lm_albania <- lm(lifeExp ~ year1950, 
                 data = by_country$data[[2]])
lm_algeria <- lm(lifeExp ~ year1950, 
                 data = by_country$data[[3]])

🙈 But, we are copying and pasting this code more than twice… Is there a better way?

A case for map???

map(<data object>, <function>)

1. Write a function to fit lm to all countries and map the function to data column of by_country

fit_lm <- function(x){
  lm(lifeExp ~ year1950, data = x) 
}

mapped_lm <- map(by_country$data, fit_lm)

head(mapped_lm)
[[1]]

Call:
lm(formula = lifeExp ~ year1950, data = x)

Coefficients:
(Intercept)     year1950  
    29.3566       0.2753  


[[2]]

Call:
lm(formula = lifeExp ~ year1950, data = x)

Coefficients:
(Intercept)     year1950  
    58.5598       0.3347  


[[3]]

Call:
lm(formula = lifeExp ~ year1950, data = x)

Coefficients:
(Intercept)     year1950  
    42.2364       0.5693  


[[4]]

Call:
lm(formula = lifeExp ~ year1950, data = x)

Coefficients:
(Intercept)     year1950  
    31.7080       0.2093  


[[5]]

Call:
lm(formula = lifeExp ~ year1950, data = x)

Coefficients:
(Intercept)     year1950  
    62.2250       0.2317  


[[6]]

Call:
lm(formula = lifeExp ~ year1950, data = x)

Coefficients:
(Intercept)     year1950  
    67.9451       0.2277  

2. Map inside the data?

country_model <- by_country %>% 
                    mutate(model = map(data, function(x){
                      lm(lifeExp ~ year1950, data = x)
                      })
                      )

country_model
# A tibble: 142 × 4
# Groups:   country, continent [142]
   country     continent data              model 
   <fct>       <fct>     <list>            <list>
 1 Afghanistan Asia      <tibble [12 × 2]> <lm>  
 2 Albania     Europe    <tibble [12 × 2]> <lm>  
 3 Algeria     Africa    <tibble [12 × 2]> <lm>  
 4 Angola      Africa    <tibble [12 × 2]> <lm>  
 5 Argentina   Americas  <tibble [12 × 2]> <lm>  
 6 Australia   Oceania   <tibble [12 × 2]> <lm>  
 7 Austria     Europe    <tibble [12 × 2]> <lm>  
 8 Bahrain     Asia      <tibble [12 × 2]> <lm>  
 9 Bangladesh  Asia      <tibble [12 × 2]> <lm>  
10 Belgium     Europe    <tibble [12 × 2]> <lm>  
# ℹ 132 more rows

3. Case for map (shorthand function)

country_model <- by_country %>% 
  mutate(model = map(data, ~lm(lifeExp ~ year1950,
                               data = .)))

country_model
# A tibble: 142 × 4
# Groups:   country, continent [142]
   country     continent data              model 
   <fct>       <fct>     <list>            <list>
 1 Afghanistan Asia      <tibble [12 × 2]> <lm>  
 2 Albania     Europe    <tibble [12 × 2]> <lm>  
 3 Algeria     Africa    <tibble [12 × 2]> <lm>  
 4 Angola      Africa    <tibble [12 × 2]> <lm>  
 5 Argentina   Americas  <tibble [12 × 2]> <lm>  
 6 Australia   Oceania   <tibble [12 × 2]> <lm>  
 7 Austria     Europe    <tibble [12 × 2]> <lm>  
 8 Bahrain     Asia      <tibble [12 × 2]> <lm>  
 9 Bangladesh  Asia      <tibble [12 × 2]> <lm>  
10 Belgium     Europe    <tibble [12 × 2]> <lm>  
# ℹ 132 more rows

4. What’s the model?

country_model$model[[1]]

Call:
lm(formula = lifeExp ~ year1950, data = .)

Coefficients:
(Intercept)     year1950  
    29.3566       0.2753  

5. How do we summarise the content in the model?

tidy(country_model$model[[1]])
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   29.4      0.699       42.0 1.40e-12
2 year1950       0.275    0.0205      13.5 9.84e- 8

6. So should we repeat it for each one?

tidy(country_model$model[[1]])
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   29.4      0.699       42.0 1.40e-12
2 year1950       0.275    0.0205      13.5 9.84e- 8
tidy(country_model$model[[2]])
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   58.6      1.13        51.7 1.79e-13
2 year1950       0.335    0.0332      10.1 1.46e- 6
tidy(country_model$model[[3]])
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   42.2      0.756       55.8 8.22e-14
2 year1950       0.569    0.0221      25.7 1.81e-10

NO!! There’s a better way. 😙

country_model %>%
  mutate(tidy = map(model, tidy))
# A tibble: 142 × 5
# Groups:   country, continent [142]
   country     continent data              model  tidy            
   <fct>       <fct>     <list>            <list> <list>          
 1 Afghanistan Asia      <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 2 Albania     Europe    <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 3 Algeria     Africa    <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 4 Angola      Africa    <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 5 Argentina   Americas  <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 6 Australia   Oceania   <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 7 Austria     Europe    <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 8 Bahrain     Asia      <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
 9 Bangladesh  Asia      <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
10 Belgium     Europe    <tibble [12 × 2]> <lm>   <tibble [2 × 5]>
# ℹ 132 more rows

Data Wrangle and tidy:

country_coefs <- country_model %>%
                    mutate(tidy = map(model, tidy)) %>%
                    unnest(tidy) %>%
                    select(country, continent, term, estimate)

country_coefs
# A tibble: 284 × 4
# Groups:   country, continent [142]
   country     continent term        estimate
   <fct>       <fct>     <chr>          <dbl>
 1 Afghanistan Asia      (Intercept)   29.4  
 2 Afghanistan Asia      year1950       0.275
 3 Albania     Europe    (Intercept)   58.6  
 4 Albania     Europe    year1950       0.335
 5 Algeria     Africa    (Intercept)   42.2  
 6 Algeria     Africa    year1950       0.569
 7 Angola      Africa    (Intercept)   31.7  
 8 Angola      Africa    year1950       0.209
 9 Argentina   Americas  (Intercept)   62.2  
10 Argentina   Americas  year1950       0.232
# ℹ 274 more rows
tidy_country_coefs <- country_coefs %>%
                          pivot_wider(id_cols = c(country, continent), 
                                      names_from =  term,
                                      values_from = estimate) %>%
                          rename(intercept = `(Intercept)`)

tidy_country_coefs
# A tibble: 142 × 4
# Groups:   country, continent [142]
   country     continent intercept year1950
   <fct>       <fct>         <dbl>    <dbl>
 1 Afghanistan Asia           29.4    0.275
 2 Albania     Europe         58.6    0.335
 3 Algeria     Africa         42.2    0.569
 4 Angola      Africa         31.7    0.209
 5 Argentina   Americas       62.2    0.232
 6 Australia   Oceania        67.9    0.228
 7 Austria     Europe         66.0    0.242
 8 Bahrain     Asia           51.8    0.468
 9 Bangladesh  Asia           35.1    0.498
10 Belgium     Europe         67.5    0.209
# ℹ 132 more rows

7. Check for Australia

tidy_country_coefs %>%
  filter(country == "Australia")
# A tibble: 1 × 4
# Groups:   country, continent [1]
  country   continent intercept year1950
  <fct>     <fct>         <dbl>    <dbl>
1 Australia Oceania        67.9    0.228

Section C: Extension Exercise

  • Fit the models to all countries.
  • Pick your favourite country (other than Australia), print the coefficients, and make a hand sketch of the the model fit.

1. Plot all the models

country_aug <- country_model %>% 
  mutate(augmented = map(model, augment)) %>%
  unnest(augmented)

country_aug
# A tibble: 1,704 × 12
# Groups:   country, continent [142]
   country     continent data     model  lifeExp year1950 .fitted  .resid   .hat
   <fct>       <fct>     <list>   <list>   <dbl>    <dbl>   <dbl>   <dbl>  <dbl>
 1 Afghanistan Asia      <tibble> <lm>      28.8        2    29.9 -1.11   0.295 
 2 Afghanistan Asia      <tibble> <lm>      30.3        7    31.3 -0.952  0.225 
 3 Afghanistan Asia      <tibble> <lm>      32.0       12    32.7 -0.664  0.169 
 4 Afghanistan Asia      <tibble> <lm>      34.0       17    34.0 -0.0172 0.127 
 5 Afghanistan Asia      <tibble> <lm>      36.1       22    35.4  0.674  0.0991
 6 Afghanistan Asia      <tibble> <lm>      38.4       27    36.8  1.65   0.0851
 7 Afghanistan Asia      <tibble> <lm>      39.9       32    38.2  1.69   0.0851
 8 Afghanistan Asia      <tibble> <lm>      40.8       37    39.5  1.28   0.0991
 9 Afghanistan Asia      <tibble> <lm>      41.7       42    40.9  0.754  0.127 
10 Afghanistan Asia      <tibble> <lm>      41.8       47    42.3 -0.534  0.169 
# ℹ 1,694 more rows
# ℹ 3 more variables: .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>
p1 <- gapminder %>% 
  ggplot(aes(year, lifeExp, group = country)) +
    geom_line(alpha = 1/3) + ggtitle("Data")

p2 <- ggplot(country_aug) + 
  geom_line(aes(x = year1950 + 1950, 
                y = .fitted, 
                group = country), 
            alpha = 1/3) +
  xlab("year") +
  ggtitle("Fitted models")
library(gridExtra)
grid.arrange(p1, p2, ncol=2)

2. Plot all the model coefficients

p <- ggplot(tidy_country_coefs, 
            aes(x = intercept, 
                y = year1950, 
                colour = continent, 
                label = country)) +
  geom_point(alpha = 0.5, 
             size = 2) +
  scale_color_brewer(palette = "Dark2")
library(plotly)
ggplotly(p)

Let’s summarise the information we learnt from the model coefficients.

  • Generally, the relationship is negative: this means that if a country starts with a high intercept, it tends to have lower rate of increase.
  • There is a difference across the continents: Countries in Europe and Oceania tend to start with a higher life expectancy and increased; countries in Asia and America tended to start lower but have high rates of improvement; Africa tends to start lower and have a huge range in rate of change.
  • Three countries had negative growth in life expectancy: Rwand, Zimbabwe, Zambia

Model diagnostics by country

country_glance <- country_model %>% 
  mutate(glance = map(model, glance)) %>%
  unnest(glance)

country_glance
# A tibble: 142 × 16
# Groups:   country, continent [142]
   country     continent data     model  r.squared adj.r.squared sigma statistic
   <fct>       <fct>     <list>   <list>     <dbl>         <dbl> <dbl>     <dbl>
 1 Afghanistan Asia      <tibble> <lm>       0.948         0.942 1.22      181. 
 2 Albania     Europe    <tibble> <lm>       0.911         0.902 1.98      102. 
 3 Algeria     Africa    <tibble> <lm>       0.985         0.984 1.32      662. 
 4 Angola      Africa    <tibble> <lm>       0.888         0.877 1.41       79.1
 5 Argentina   Americas  <tibble> <lm>       0.996         0.995 0.292    2246. 
 6 Australia   Oceania   <tibble> <lm>       0.980         0.978 0.621     481. 
 7 Austria     Europe    <tibble> <lm>       0.992         0.991 0.407    1261. 
 8 Bahrain     Asia      <tibble> <lm>       0.967         0.963 1.64      291. 
 9 Bangladesh  Asia      <tibble> <lm>       0.989         0.988 0.977     930. 
10 Belgium     Europe    <tibble> <lm>       0.995         0.994 0.293    1822. 
# ℹ 132 more rows
# ℹ 8 more variables: p.value <dbl>, df <dbl>, logLik <dbl>, AIC <dbl>,
#   BIC <dbl>, deviance <dbl>, df.residual <int>, nobs <int>

3. Plot the \(R^2\) values using a histogram.

ggplot(country_glance, 
       aes(x = r.squared)) + 
  geom_histogram()

4. Countries with worst fit

Examine the countries with the worst fit, countries with \(R^2<0.45\), by making scatterplots of the data, with the linear model overlaid.

badfit <- country_glance %>% filter(r.squared <= 0.45)

gap_bad <- gap %>% filter(country %in% badfit$country)

gg_bad_fit <-
ggplot(data = gap_bad, 
       aes(x = year, 
           y = lifeExp)) + 
         geom_point() +
  facet_wrap(~country) +
  scale_x_continuous(breaks = seq(1950,2000,10), 
                     labels = c("1950", "60","70", "80","90","2000")) +
  geom_smooth(method = "lm", 
              se = FALSE)
gg_bad_fit

Each of these countries was moving on a nice trajectory of increasing life expectancy, but later suffered a big dip during the time period.